Turbine Collision Risk: Exploration Circular Distributions and Acoustic Monitoring Area

Cédric Scherer https://cedricscherer.com (Leibniz Institute for Zoo and Wildlife Research)https://www.izw-berlin.de/en/home.html , Christian Voigt (Leibniz Institute for Zoo and Wildlife Research)https://www.izw-berlin.de/en/home.html
2021-11-25

Preparation

Show code
library(tidyverse)
library(ggtext)
library(spatstat.core)
library(patchwork)

theme_set(theme_light(base_size = 10, base_family = "Open Sans"))
theme_update(
  panel.grid.major = element_line(size = .3),
  panel.grid.minor = element_blank(),
  plot.title = element_markdown(size = 18),
  plot.title.position = "plot",
  plot.margin = margin(rep(25, 4))
)

Research Question

Is the predictive power of the current acoustic monitoring sufficient enough to estimate the true number of killed bats?

Predictions

The predictive power of the acoustic monitoring for the extrapolation of the expected number of stroke victims is decreasing

Open Questions

Inputs

We have the following parameters we can vary:

Show code
## radius circle (~ length rotor blades)
radius_rotor <- 1 

## % area covered by acoustic monitoring
proportion_covered <- .05 

## radius acoustic monitoring via % area covered + length rotor blades
area_rotor <- pi * radius_rotor^2
area_monitored <- area_rotor * proportion_covered
radius_monitored <- sqrt(area_monitored / pi)

## mask for window
circle <- disc(
  radius = radius_rotor, 
  centre = c(0, 0), 
  mask = FALSE, 
  npoly = 5000
)

## number of points (~ bats)
n <- 500

Circular Distributions

Generate Point Patterns

Show code
## source sunflower funciton for uniform distributions
source(here::here("R", "src", "util-uniform-distribution.R"))

## fix random number genreration to ensure reproducibility
set.seed(12345)

## random
pp_rand <- spatstat.core::rpoint(n, win = circle)
df_rand <- as.data.frame(pp_rand) %>% mutate(dist = x^2 + y^2)


## uniform
pp_unif <- spatstat.core::runifpoint(n, win = circle)
df_unif <- as.data.frame(pp_unif) %>% mutate(dist = x^2 + y^2)

df_sunf <- sunflower(n, 1, 'planar') %>% mutate(dist = x^2 + y^2)


## directional: inner > outer
pp_inn <- rpoint(n, function(x,y) {1 - (abs(-x^2 - y^2)) + .01}, win = circle)
df_inn <- as.data.frame(pp_inn) %>% mutate(dist = x^2 + y^2)


## directional: outer > inner
pp_out <- rpoint(n, function(x,y) {(abs(x^2 + y^2)) + .01}, win = circle)
df_out <- as.data.frame(pp_out) %>% mutate(dist = x^2 + y^2)


## directional: bottom > top
pp_btm1 <- rpoint(n, function(x,y) {100 * exp(-1*y)}, win = circle)
df_btm1 <- as.data.frame(pp_btm1) %>% mutate(dist = x^2 + y^2)

pp_btm2 <- rpoint(n, function(x,y) {100 * exp(-2*y)}, win = circle)
df_btm2 <- as.data.frame(pp_btm2) %>% mutate(dist = x^2 + y^2)

pp_btm3 <- rpoint(n, function(x,y) {100 * exp(-3*y)}, win = circle)
df_btm3 <- as.data.frame(pp_btm3) %>% mutate(dist = x^2 + y^2)

pp_btm4 <- rpoint(n, function(x,y) {100 * exp(-4*y)}, win = circle)
df_btm4 <- as.data.frame(pp_btm4) %>% mutate(dist = x^2 + y^2)

pp_btm5 <- rpoint(n, function(x,y) {100 * exp(-5*y)}, win = circle)
df_btm5 <- as.data.frame(pp_btm5) %>% mutate(dist = x^2 + y^2)

pp_btm6 <- rpoint(n, function(x,y) {100 * exp(-6*y)}, win = circle)
df_btm6 <- as.data.frame(pp_btm6) %>% mutate(dist = x^2 + y^2)

pp_btm7 <- rpoint(n, function(x,y) {100 * exp(-7*y)}, win = circle)
df_btm7 <- as.data.frame(pp_btm7) %>% mutate(dist = x^2 + y^2)

pp_btm8 <- rpoint(n, function(x,y) {100 * exp(-8*y)}, win = circle)
df_btm8 <- as.data.frame(pp_btm8) %>% mutate(dist = x^2 + y^2)

pp_btm9 <- rpoint(n, function(x,y) {100 * exp(-9*y)}, win = circle)
df_btm9 <- as.data.frame(pp_btm9) %>% mutate(dist = x^2 + y^2)

pp_btm10 <- rpoint(n, function(x,y) {100 * exp(-10*y)}, win = circle)
df_btm10 <- as.data.frame(pp_btm10) %>% mutate(dist = x^2 + y^2)

Plots Point Patterns

Show code
plot_dist <- function(data, title, color = TRUE) {

  base <- 
    ggplot(as.data.frame(circle), aes(x, y)) + 
    geom_path(color = "black", alpha = .33, size = 2) + 
    coord_fixed() + 
    scale_x_continuous(limits = c(-1, 1)) +
    scale_y_continuous(limits = c(-1, 1))
  
  if (isTRUE(color)) {
     base + 
      geom_point(data = data, aes(color = dist), alpha = .7) +
      scale_color_gradient(low = "#003F4C", high = "#B200B2", guide = "none") + 
      ggtitle(title)
  } else {
    base + 
      geom_point(data = data, alpha = .7) + 
      ggtitle(title)
  }
}

Distribution Scenarios

Show code
## no color versions ###########################################################

## random
gg_rand_nc <- plot_dist(df_rand, "**Random**", color = FALSE)

## uniform
#gg_unif_nc <- plot_dist(df_unif, "**Uniform** (?)", color = FALSE)
gg_sunf_nc <- plot_dist(df_sunf, "**Uniform** (sunflower)", color = FALSE)  

## directional: inner > outer
gg_inn_nc <- plot_dist(df_inn, "**Inner to Outer**", color = FALSE)

## directional: outer > inner
gg_out_nc <- plot_dist(df_out, "**Outer to Inner**", color = FALSE)               

## directional: bottom > top
gg_btm1_nc  <- plot_dist(df_btm1,  "**Bottom to Top** (exp(-y))", color = FALSE)
gg_btm2_nc  <- plot_dist(df_btm2,  "**Bottom to Top** (exp(-2*y))", color = FALSE)
gg_btm3_nc  <- plot_dist(df_btm3,  "**Bottom to Top** (exp(-3*y))", color = FALSE)
gg_btm4_nc  <- plot_dist(df_btm4,  "**Bottom to Top** (exp(-4*y))", color = FALSE)
gg_btm5_nc  <- plot_dist(df_btm5,  "**Bottom to Top** (exp(-5*y))", color = FALSE)
gg_btm6_nc  <- plot_dist(df_btm6,  "**Bottom to Top** (exp(-6*y))", color = FALSE)
gg_btm7_nc  <- plot_dist(df_btm7,  "**Bottom to Top** (exp(-7*y))", color = FALSE)
gg_btm8_nc  <- plot_dist(df_btm8,  "**Bottom to Top** (exp(-8*y))", color = FALSE)
gg_btm9_nc  <- plot_dist(df_btm9,  "**Bottom to Top** (exp(-9*y))", color = FALSE)
gg_btm10_nc <- plot_dist(df_btm10, "**Bottom to Top** (exp(-10*y))", color = FALSE)


## panel
(gg_sunf_nc | gg_btm2_nc | gg_inn_nc) / (gg_rand_nc | gg_btm5_nc | gg_out_nc)
Show code
ggsave("plots/circular_distributions_nc.png", width = 15, height = 10, dpi = 700)

Distribution Scenarios Colored by Distance

Show code
## colored versions ############################################################

## random
gg_rand <- plot_dist(df_rand, "**Random**")

## uniform
#gg_unif <- plot_dist(df_unif, "**Uniform** (?)")
gg_sunf <- plot_dist(df_sunf, "**Uniform** (sunflower)")  

## directional: inner > outer
gg_inn <- plot_dist(df_inn, "**Inner to Outer**")

## directional: outer > inner
gg_out <- plot_dist(df_out, "**Outer to Inner**")

## directional: bottom > top
gg_btm1  <- plot_dist(df_btm1,  "**Bottom to Top** (exp(-y))")
gg_btm2  <- plot_dist(df_btm2,  "**Bottom to Top** (exp(-2*y))")
gg_btm3  <- plot_dist(df_btm3,  "**Bottom to Top** (exp(-3*y))")
gg_btm4  <- plot_dist(df_btm4,  "**Bottom to Top** (exp(-4*y))")
gg_btm5  <- plot_dist(df_btm5,  "**Bottom to Top** (exp(-5*y))")
gg_btm6  <- plot_dist(df_btm6,  "**Bottom to Top** (exp(-6*y))")
gg_btm7  <- plot_dist(df_btm7,  "**Bottom to Top** (exp(-7*y))")
gg_btm8  <- plot_dist(df_btm8,  "**Bottom to Top** (exp(-8*y))")
gg_btm9  <- plot_dist(df_btm9,  "**Bottom to Top** (exp(-9*y))")
gg_btm10 <- plot_dist(df_btm10, "**Bottom to Top** (exp(-10*y))")


## panel
(gg_sunf | gg_btm2 | gg_inn) / (gg_rand | gg_btm5 | gg_out)
Show code
ggsave("plots/circular_distributions.png", width = 15, height = 10, dpi = 700)

Variations of Bottom to Top Distributions

Show code
## panel bottom > top
(gg_btm1_nc | gg_btm2_nc | gg_btm3_nc | gg_btm4_nc | gg_btm5_nc) / (gg_btm6_nc | gg_btm7_nc | gg_btm8_nc | gg_btm9_nc | gg_btm10_nc)
Show code
ggsave("plots/circular_distributions_btms_nc.png", width = 24, height = 10, dpi = 700)

Distribution Scenarios Colored by Distance

Show code
## panel bottom > top
(gg_btm1 | gg_btm2 | gg_btm3 | gg_btm4 | gg_btm5) / (gg_btm6 | gg_btm7 | gg_btm8 | gg_btm9 | gg_btm10)
Show code
ggsave("plots/circular_distributions_btms.png", width = 24, height = 10, dpi = 700)

Plots Proportion Monitored

Show code
plot_monitoring <- function(gg, data, circle_monitored) {
  gg + 
    geom_path(data = as.data.frame(circle_monitored), color = "orange2", alpha = .2, size = 1.5) +
    geom_point(data = filter(data, dist < radius_monitored^2), color = "orange2")
}
Show code
plot_monitoring_all <- function(proportion_covered, df_rand, df_sunf, df_inn, df_out, df_btm_2, df_btm_5) {
  
  ## radius acoustic monitoring via % area covered + length rotor blades
  area_rotor <- pi * radius_rotor^2
  area_monitored <- area_rotor * proportion_covered * 2   
  ## times 2 because we only use the lower half of the circle = half of the area later
  radius_monitored <- sqrt(area_monitored / pi)
  
  circle_monitored <- disc(
    radius = radius_monitored,
    centre = c(0, 0), 
    mask = FALSE, 
    npoly = 5000
  )
  
  circle_half_monited <-
    as.data.frame(circle_monitored) %>% 
    filter(y <= 0)
    
  ## random
  gg_rand_r <- gg_rand_nc + 
    geom_path(data = circle_half_monited, color = "orange2", alpha = .2, size = 1.5) +
    geom_point(data = filter(df_rand, dist < radius_monitored^2 & y <= 0), color = "orange2")
  
  ## uniform
  # gg_unif_r <- gg_unif_nc + 
  #   geom_path(data = as.data.frame(circle_monitored), color = "orange2", alpha = .2, size = 1.5) +
  #   geom_point(data = filter(df_unif, dist < radius_monitored^2), color = "orange2")
  
  gg_sunf_r <- gg_sunf_nc + 
    geom_path(data = as.data.frame(circle_monitored), color = "orange2", alpha = .2, size = 1.5) +
    geom_point(data = filter(df_sunf, dist < radius_monitored^2), color = "orange2")

  ## directional: inner > outer
  gg_inn_r <- gg_inn_nc + 
    geom_path(data = as.data.frame(circle_monitored), color = "orange2", alpha = .2, size = 1.5) +
    geom_point(data = filter(df_inn, dist < radius_monitored^2), color = "orange2")
  
  ## directional: outer > inner
  gg_out_r <- gg_out_nc + 
    geom_path(data = as.data.frame(circle_monitored), color = "orange2", alpha = .2, size = 1.5) +
    geom_point(data = filter(df_out, dist < radius_monitored^2), color = "orange2")
  
  
  ## bottom > top
  gg_btm2_r <- gg_btm2_nc + 
    geom_path(data = as.data.frame(circle_monitored), color = "orange2", alpha = .2, size = 1.5) +
    geom_point(data = filter(df_btm2, dist < radius_monitored^2), color = "orange2")
  
  gg_btm5_r <- gg_btm5_nc + 
    geom_path(data = as.data.frame(circle_monitored), color = "orange2", alpha = .2, size = 1.5) +
    geom_point(data = filter(df_btm5, dist < radius_monitored^2), color = "orange2")
  
  
  ## panel
  p <- (gg_sunf_r | gg_btm2_r | gg_inn_r) / (gg_rand_r | gg_btm5_r | gg_out_r)
}

Monitoring Covers 50%

Show code
(panel_monitored <- plot_monitoring_all(
  proportion_covered = .5, df_rand, df_sunf, df_inn, df_out, df_btm_2, df_btm_5
))
Show code
ggsave("plots/circular_distributions_area_monitored_50.png", width = 15, height = 10, dpi = 700)

Monitoring Covers 40%

Show code
(panel_monitored <- plot_monitoring_all(
  proportion_covered = .4, df_rand, df_sunf, df_inn, df_out, df_btm_2, df_btm_5
))
Show code
ggsave("plots/circular_distributions_area_monitored_40.png", width = 15, height = 10, dpi = 700)

Monitoring Covers 30%

Show code
(panel_monitored <- plot_monitoring_all(
  proportion_covered = .3, df_rand, df_sunf, df_inn, df_out, df_btm_2, df_btm_5
))
Show code
ggsave("plots/circular_distributions_area_monitored_30.png", width = 15, height = 10, dpi = 700)

Monitoring Covers 20%

Show code
(panel_monitored <- plot_monitoring_all(
  proportion_covered = .2, df_rand, df_sunf, df_inn, df_out, df_btm_2, df_btm_5
))
Show code
ggsave("plots/circular_distributions_area_monitored_20.png", width = 15, height = 10, dpi = 700)

Monitoring Covers 10%

Show code
(panel_monitored <- plot_monitoring_all(
  proportion_covered = .1, df_rand, df_sunf, df_inn, df_out, df_btm_2, df_btm_5
))
Show code
ggsave("plots/circular_distributions_area_monitored_10.png", width = 15, height = 10, dpi = 700)

Monitoring Covers 5%

Show code
(panel_monitored <- plot_monitoring_all(
  proportion_covered = .05, df_rand, df_sunf, df_inn, df_out, df_btm_2, df_btm_5
))
Show code
ggsave("plots/circular_distributions_area_monitored_05.png", width = 15, height = 10, dpi = 700)

Session Info
Show code
[1] "2021-11-25 17:22:37 CET"
Show code
R version 4.1.0 (2021-05-18)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19043)

Matrix products: default

locale:
[1] LC_COLLATE=German_Germany.1252  LC_CTYPE=German_Germany.1252   
[3] LC_MONETARY=German_Germany.1252 LC_NUMERIC=C                   
[5] LC_TIME=German_Germany.1252    
system code page: 65001

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods  
[7] base     

other attached packages:
 [1] patchwork_1.1.1     spatstat.core_2.3-0 rpart_4.1-15       
 [4] nlme_3.1-152        spatstat.geom_2.3-0 spatstat.data_2.1-0
 [7] ggtext_0.1.1        forcats_0.5.1       stringr_1.4.0      
[10] dplyr_1.0.7         purrr_0.3.4         readr_2.0.2        
[13] tidyr_1.1.4         tibble_3.1.5        ggplot2_3.3.5      
[16] tidyverse_1.3.1    

loaded via a namespace (and not attached):
 [1] fs_1.5.0              spatstat.sparse_2.0-0 lubridate_1.8.0      
 [4] httr_1.4.2            rprojroot_2.0.2       tools_4.1.0          
 [7] backports_1.2.1       bslib_0.3.1           utf8_1.2.2           
[10] R6_2.5.1              DBI_1.1.1             mgcv_1.8-35          
[13] colorspace_2.0-2      withr_2.4.2           tidyselect_1.1.1     
[16] downlit_0.2.1         compiler_4.1.0        textshaping_0.3.6    
[19] cli_3.0.0             rvest_1.0.2           xml2_1.3.2           
[22] labeling_0.4.2        sass_0.4.0            scales_1.1.1         
[25] goftest_1.2-3         systemfonts_1.0.3     digest_0.6.28        
[28] spatstat.utils_2.2-0  rmarkdown_2.11        pkgconfig_2.0.3      
[31] htmltools_0.5.2       highr_0.9             dbplyr_2.1.1         
[34] fastmap_1.1.0         rlang_0.4.12          readxl_1.3.1         
[37] rstudioapi_0.13       farver_2.1.0          jquerylib_0.1.4      
[40] generics_0.1.0        jsonlite_1.7.2        distill_1.3          
[43] magrittr_2.0.1        Matrix_1.3-3          Rcpp_1.0.7           
[46] munsell_0.5.0         fansi_0.5.0           abind_1.4-5          
[49] lifecycle_1.0.1       stringi_1.7.5         yaml_2.2.1           
[52] grid_4.1.0            crayon_1.4.1          deldir_1.0-5         
[55] lattice_0.20-44       haven_2.4.3           splines_4.1.0        
[58] tensor_1.5            gridtext_0.1.4        hms_1.1.1            
[61] knitr_1.36            pillar_1.6.4          markdown_1.1         
[64] reprex_2.0.1          glue_1.4.2            evaluate_0.14        
[67] modelr_0.1.8          vctrs_0.3.8           tzdb_0.1.2           
[70] cellranger_1.1.0      gtable_0.3.0          polyclip_1.10-0      
[73] assertthat_0.2.1      xfun_0.27             broom_0.7.9          
[76] ragg_1.1.3            ellipsis_0.3.2        here_1.0.1